pacman::p_load(dplyr,tidyr, rcompanion, readxl,ggplot2, ggstatsplot,brms, bayesplot, rstan, tidybayes,ggdist)
df<-read_excel("Food_Habits_109_sex_answers.xlsx")%>%
dplyr::select(`Group (as in the html file: 1, CoV_smell_taste_disturbances; 2, CoV_no_smell_taste_disturbances; 3, healthCtr)`,Number, Sex,
`How much salt do you use in your meals?`,
`How much bitter additives do you use in your meals?`,
`How much sugar/sweeteners do you use in your meals?`,
`How much sour additives (e.g. citric acid, vinegar) do you use in your meals?`)%>%
rename("Group"="Group (as in the html file: 1, CoV_smell_taste_disturbances; 2, CoV_no_smell_taste_disturbances; 3, healthCtr)",
"ID"="Number",
"salt"="How much salt do you use in your meals?",
"bitter"="How much bitter additives do you use in your meals?",
"sweet"="How much sugar/sweeteners do you use in your meals?",
"sour"="How much sour additives (e.g. citric acid, vinegar) do you use in your meals?")%>%
mutate(Group=recode(Group, "1"="CoV_smell_taste_disturbances",
"2"="CoV_no_smell_taste_disturbances","3"="healthCtr"))%>%
drop_na(ID)%>%
mutate_all(~ifelse(.=="None",0, ifelse(.=="Little",1,ifelse(.=="VeryLittle", 2,
ifelse(.=="Moderately",3, ifelse(.=="Much",4, ifelse(.=="VeryMuch", 5,.)))))))%>%
mutate(salt=as.numeric(salt),
sweet=as.numeric(sweet),
bitter=as.numeric(bitter),
sour=as.numeric(sour))%>%
glimpse()
## Rows: 71
## Columns: 7
## $ Group <chr> "CoV_smell_taste_disturbances", "CoV_smell_taste_disturbances",…
## $ ID <chr> "GDK003", "GDK004", "GDK007", "GDK011", "GDK012", "GDK021", "GD…
## $ Sex <chr> "M", "M", "F", "F", "F", "F", "F", "M", "F", "F", "F", "F", "M"…
## $ salt <dbl> 3, 3, 3, 3, 3, 3, 4, 4, 2, 4, 2, 3, 4, 4, 3, 1, 3, 3, 3, 3, 4, …
## $ bitter <dbl> 3, 1, 4, 3, 3, 5, 4, 1, 1, 4, 3, 0, 3, 3, 3, 3, 2, 1, 1, 1, 4, …
## $ sweet <dbl> 3, 5, 3, 2, 1, 1, 4, 3, 2, 2, 3, 0, 2, 2, 3, 1, 3, 0, 3, 4, 1, …
## $ sour <dbl> 3, 3, 1, 1, 2, 0, 3, 3, 3, 3, 1, 2, 3, 3, 1, 2, 2, 2, 1, 2, 3, …
summary(df)
## Group ID Sex salt
## Length:71 Length:71 Length:71 Min. :1.000
## Class :character Class :character Class :character 1st Qu.:2.000
## Mode :character Mode :character Mode :character Median :3.000
## Mean :2.662
## 3rd Qu.:3.000
## Max. :4.000
## bitter sweet sour
## Min. :0.000 Min. :0.000 Min. :0.00
## 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:1.00
## Median :3.000 Median :2.000 Median :2.00
## Mean :2.507 Mean :2.324 Mean :2.07
## 3rd Qu.:3.000 3rd Qu.:3.000 3rd Qu.:3.00
## Max. :5.000 Max. :5.000 Max. :4.00
ggbetweenstats(
data = df,
x = Group,
y = salt,
title="Salt"
)
ggbetweenstats(
data = df,
x = Group,
y = sweet,
title="Sweet"
)
ggbetweenstats(
data = df,
x = Group,
y = bitter,
title="bitter"
)
ggbetweenstats(
data = df,
x = Group,
y = sour,
title="Sour"
)
# subset data
df1<-subset(df, Group !="CoV_no_smell_taste_disturbances")
ggbetweenstats(
data = df1,
x = Group,
y = salt,
title="Salt : CoV_smell_taste_disturbances vs healthCtr"
)
ggbetweenstats(
data = df1,
x = Group,
y = sweet,
title="Sweet : CoV_smell_taste_disturbances vs healthCtr"
)
ggbetweenstats(
data = df1,
x = Group,
y = bitter,
title="bitter : CoV_smell_taste_disturbances vs healthCtr"
)
ggbetweenstats(
data = df1,
x = Group,
y = sour,
title="Sour : CoV_smell_taste_disturbances vs healthCtr"
)
#salt
scheirerRayHare(salt~Group+Sex, data=df)
##
## DV: salt
## Observations: 71
## D: 0.8826626
## MS total: 426
## Df Sum Sq H p.value
## Group 2 369.9 0.9836 0.61151
## Sex 1 32.9 0.0875 0.76737
## Group:Sex 2 1930.0 5.1328 0.07681
## Residuals 65 23988.2
scheirerRayHare(salt~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
##
## DV: salt
## Observations: 64
## D: 0.8804716
## MS total: 346.6667
## Df Sum Sq H p.value
## Group 1 263.4 0.8629 0.35293
## Sex 1 0.3 0.0009 0.97566
## Group:Sex 1 1381.3 4.5255 0.03339
## Residuals 60 17583.2
summary(aov(salt~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 0.67 0.673 0.733 0.395
## Sex 1 0.00 0.003 0.003 0.956
## Group:Sex 1 4.33 4.325 4.709 0.034 *
## Residuals 60 55.11 0.918
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sweet
scheirerRayHare(sweet~Group+Sex, data=df)
##
## DV: sweet
## Observations: 71
## D: 0.9414655
## MS total: 426
## Df Sum Sq H p.value
## Group 2 341.5 0.8515 0.65328
## Sex 1 239.7 0.5976 0.43948
## Group:Sex 2 1500.1 3.7403 0.15410
## Residuals 65 26005.3
scheirerRayHare(sweet~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
##
## DV: sweet
## Observations: 64
## D: 0.9405678
## MS total: 346.6667
## Df Sum Sq H p.value
## Group 1 260.6 0.7991 0.37136
## Sex 1 80.7 0.2474 0.61888
## Group:Sex 1 1066.2 3.2698 0.07057
## Residuals 60 19146.9
summary(aov(sweet~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 0.95 0.954 0.621 0.434
## Sex 1 0.69 0.691 0.450 0.505
## Group:Sex 1 5.93 5.932 3.862 0.054 .
## Residuals 60 92.17 1.536
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#bitter
scheirerRayHare(bitter~Group+Sex, data=df)
##
## DV: bitter
## Observations: 71
## D: 0.9156439
## MS total: 426
## Df Sum Sq H p.value
## Group 2 981.6 2.51663 0.28413
## Sex 1 16.5 0.04222 0.83719
## Group:Sex 2 553.4 1.41882 0.49194
## Residuals 65 25707.5
scheirerRayHare(bitter~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
##
## DV: bitter
## Observations: 64
## D: 0.9123397
## MS total: 346.6667
## Df Sum Sq H p.value
## Group 1 0.3 0.00083 0.97697
## Sex 1 33.3 0.10528 0.74559
## Group:Sex 1 418.6 1.32341 0.24998
## Residuals 60 19473.6
summary(aov(bitter~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 0.01 0.0135 0.009 0.926
## Sex 1 0.18 0.1831 0.120 0.731
## Group:Sex 1 1.56 1.5638 1.022 0.316
## Residuals 60 91.85 1.5308
#sour
scheirerRayHare(sour~Group+Sex, data=df)
##
## DV: sour
## Observations: 71
## D: 0.9074111
## MS total: 426
## Df Sum Sq H p.value
## Group 2 439.0 1.1356 0.56678
## Sex 1 1273.9 3.2955 0.06947
## Group:Sex 2 1031.2 2.6677 0.26346
## Residuals 65 24141.3
scheirerRayHare(sour~Group+Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances"))
##
## DV: sour
## Observations: 64
## D: 0.8957875
## MS total: 346.6667
## Df Sum Sq H p.value
## Group 1 269.3 0.86731 0.35170
## Sex 1 495.4 1.59527 0.20658
## Group:Sex 1 480.8 1.54812 0.21341
## Residuals 60 18284.3
summary(aov(sour~Group*Sex, data=subset(df, Group != "CoV_no_smell_taste_disturbances")))
## Df Sum Sq Mean Sq F value Pr(>F)
## Group 1 1.17 1.168 1.067 0.306
## Sex 1 1.66 1.665 1.521 0.222
## Group:Sex 1 1.42 1.417 1.294 0.260
## Residuals 60 65.69 1.095
Among several advanced statistical methods and resampling
techniques, Bayesian can provide more robust insights. We use Bayesian
model to perform analyses of the two-way (Group * Sex) data
# salt: almost no Group difference
model<-brm(salt~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: salt ~ Group * Sex
## Data: df1 (Number of observations: 64)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.68 0.22 2.25 3.11 1.00 3642 3130
## GrouphealthCtr 0.00 0.27 -0.55 0.55 1.00 3529 2437
## SexM 0.82 0.53 -0.23 1.87 1.00 2251 2394
## GrouphealthCtr:SexM -1.49 0.72 -2.90 -0.04 1.00 2621 2529
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.98 0.09 0.81 1.18 1.00 4036 2877
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
ggtitle("Salt: Posterior Distribution of Group healthCtr Effect")
# Posterior predictive checks
pp_check(model)+
ggtitle("Salt: Posterior Predictive Check for Group Differences")
# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)
plot(mcmc_intervals(as_draws_df(model)))+
ggtitle("Salt: 95% Credible Interval for Group healthCtr Effect")
# Difference in Group Means Plot
#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr
#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean
#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
ggtitle("Salt: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")
# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)
# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)
# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
stat_pointinterval() +
ggtitle("Salt: Predicted Group Differences by Sex")+theme_ggdist()
## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = salt, fill = Group)) +
stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
geom_jitter(width = .1, alpha = 0.3) +
ggtitle("Salt: Raincloud Plot of Group Differences")+theme_ggdist()+
theme(legend.position = "none")
# sweet:COVID patients reduced sugar intake
model<-brm(sweet~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: sweet ~ Group * Sex
## Data: df1 (Number of observations: 64)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.95 0.27 1.41 2.48 1.00 3275 3120
## GrouphealthCtr 0.52 0.35 -0.15 1.20 1.00 3382 2630
## SexM 1.28 0.67 -0.00 2.59 1.00 2399 3043
## GrouphealthCtr:SexM -1.75 0.90 -3.48 0.03 1.00 2263 3001
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.26 0.11 1.06 1.51 1.00 3222 2973
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
ggtitle("Sweet: Posterior Distribution of Group healthCtr Effect")
# Posterior predictive checks
pp_check(model)+
ggtitle("Sweet: Posterior Predictive Check for Group Differences")
# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)
plot(mcmc_intervals(as_draws_df(model)))+
ggtitle("Sweet: 95% Credible Interval for Group healthCtr Effect")
# Difference in Group Means Plot
#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr
#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean
#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
ggtitle("Sweet: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")
# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)
# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)
# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
stat_pointinterval() +
ggtitle("Sweet: Predicted Group Differences by Sex")+theme_ggdist()
## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = sweet, fill = Group)) +
stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
geom_jitter(width = .1, alpha = 0.3) +
ggtitle("Sweet: Raincloud Plot of Group Differences")+theme_ggdist()+
theme(legend.position = "none")
## bitter: COVID patients increase bitter additive intake!
model<-brm(bitter~Group*Sex, data=df1, family=gaussian())
# Posterior Distributions
summary(model)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: bitter ~ Group * Sex
## Data: df1 (Number of observations: 64)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.67 0.27 2.14 3.20 1.00 3203 2876
## GrouphealthCtr -0.11 0.35 -0.79 0.56 1.00 3069 2968
## SexM -0.66 0.70 -2.01 0.71 1.00 2411 2785
## GrouphealthCtr:SexM 0.90 0.91 -0.89 2.66 1.00 2277 2544
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.26 0.12 1.06 1.52 1.00 4197 2520
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Extract posterior samples
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
ggtitle("Bitter: Posterior Distribution of Group healthCtr Effect")
# Posterior predictive checks
pp_check(model)+
ggtitle("Bitter: Posterior Predictive Check for Group Differences")
# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)
plot(mcmc_intervals(as_draws_df(model)))+
ggtitle("Bitter: 95% Credible Interval for Group healthCtr Effect")
# Difference in Group Means Plot
#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr
#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean
#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
ggtitle("Bitter: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")
# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)
# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)
# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
stat_pointinterval() +
ggtitle("Bitter: Predicted Group Differences by Sex")+theme_ggdist()
## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = bitter, fill = Group)) +
stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
geom_jitter(width = .1, alpha = 0.3) +
ggtitle("Bitter: Raincloud Plot of Group Differences")+theme_ggdist()+
theme(legend.position = "none")
## sour COVID patients increase sour additive intake
model<-brm(sour~Group*Sex, data=df1, family=gaussian())
summary(model)
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: sour ~ Group * Sex
## Data: df1 (Number of observations: 64)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.05 0.23 1.58 2.50 1.00 2990 2498
## GrouphealthCtr -0.14 0.30 -0.73 0.44 1.00 2879 2438
## SexM 0.94 0.58 -0.17 2.07 1.00 2299 2331
## GrouphealthCtr:SexM -0.83 0.77 -2.29 0.63 1.00 2047 2605
##
## Further Distributional Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 1.07 0.10 0.90 1.28 1.00 4402 3014
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
post_samples <- posterior_samples(model)
# Plot posterior distributions of group effect
mcmc_dens(post_samples, pars = "b_GrouphealthCtr")+
ggtitle("Sour: Posterior Distribution of Group healthCtr Effect")
# Posterior predictive checks
pp_check(model)+
ggtitle("Sour: Posterior Predictive Check for Group Differences")
# Credible Intervals Plot
draws <- as_draws_df(model)
mcmc_trace(draws)
plot(mcmc_intervals(as_draws_df(model)))+
ggtitle("Sour: 95% Credible Interval for Group healthCtr Effect")
# Difference in Group Means Plot
#extracting group-level means from the posterior samples
CoV_smell_taste_disturbances_mean <- post_samples$b_Intercept
healthCtr_mean <- post_samples$b_Intercept + post_samples$b_GrouphealthCtr
#Difference in means
group_difference <- healthCtr_mean - CoV_smell_taste_disturbances_mean
#Plot the difference in means
mcmc_areas(as.data.frame(group_difference)) +
ggtitle("Sour: Posterior Distribution of Group Differences\n (Group_healthCtr - Group_CoV_smell_taste_disturbances)")
# Interaction Plot (for Group * Sex Interaction)
new_data <- data.frame(
Group = factor(c("CoV_smell_taste_disturbances", "healthCtr")),
Sex = factor(c("M", "F"))
)
fitted_values <- posterior_epred(model, newdata = new_data)
# Convert to a tidy format
fitted_values_df <- add_epred_draws(new_data, model)
# Plot interaction between Group and Sex
ggplot(fitted_values_df, aes(x = Group, y = .epred, color = Sex)) +
stat_pointinterval() +
ggtitle("Sour: Predicted Group Differences by Sex")+theme_ggdist()
## Raincloud Plot for Group Difference
ggplot(df1, aes(x = Group, y = sour, fill = Group)) +
stat_halfeye(adjust = .5, justification = -.2, .width = 0, point_colour = NA) +
geom_boxplot(width = .12, outlier.color = NA, alpha = 0.5) +
geom_jitter(width = .1, alpha = 0.3) +
ggtitle("Sour: Raincloud Plot of Group Differences")+theme_ggdist()+
theme(legend.position = "none")
